Direct Transition to Spatiotemporal Chaos in Low Prandtl Number Fluids 
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We present a large scale numerical simulation of three-dimensional Rayleigh-Benard convection 
near onset, under free-free boundary conditions for a fluid of Prandtl number a = 0.5. We find 
that a spatiotemporally chaotic state emerges immediately above onset, which we investigate as a 
function of the reduced control parameter. We conclude that the transition from conduction to 
spatiotemporal chaos is second order and of "mean field" character. We also present a simple theory 
for the time-averaged convective current. Finally, we show that the time-averaged structure factor 
satisfies a scaling behavior with respect to the correlation length near onset. 
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Pattern formation in non-equilibrium systems has be- 
come a major frontier area in science . The richness 
of this field has been significantly enhanced by the exis- 
tence of spatiotemporal chaos (STC) in various systems 
10-0. STC is characterized by its extensive, irregular 
dynamics in both space and time. It has been recog- 
nized in experiments [^-Q and numerical studies ||l|,||j6| 
that a large aspect ratio is essential for the occurrence of 
STC. Owing to the generic complication of its dynamics, 
theoretical understanding of STC relies heavily on some 
much simplified, mathematical models of the real sys- 
tems 0. Although much progress has been made, some 
fundamental concepts remain to be developed. 

A paradigm of pattern formation is Rayleigh-Benard 
convection (RBC) which occurs when a thin hori- 
zontal fluid layer is heated from below. In general, the 
dynamics of RBC depends on the Rayleigh number R, 
the Prandtl number a of the fluid, and the aspect ra- 
tio (size/thickness) F of the system. Busse and his col- 
laborators 1^ have studied extensively the stability do- 
main of parallel rolls as a function of wavenumber k and 
Rayleigh number R, for various a. It is well known 
that in a laterally infinite system with rigid-rigid bound- 
aries, there exists a stable, time independent parallel roll 
state near the onset of convection for all a. In the case 
of free-free boundaries at sufficiently low Prandtl num- 
bers (ct <0.543), Siggia and Zippelius |^ and Busse and 
Bolton found surprisingly that parallel rolls are un- 
stable with respect to the skewed-varicose instability im- 
mediately above onset. Busse et al. fll|] further stud- 
ied the possibility of a direction transition from conduc- 
tion to STC, but their aspect ratio (F = 8) is not large 
enough for a conclusive result. Although the free- free 
boundaries are very difficult to control for detailed ex- 
perimental studies, one experiment with such boundary 
conditions has been reported Q . 

In this paper, we present the results of a large scale 
(F = 60) numerical simulation of the three dimensional 
hydrodynamic equations, using the Boussinesq approxi- 
mation, for a low Prandtl number fluid (tT=0.5) with free- 
free boundary conditions. The same problem has been 



investigated before |T|,|l^, but the aspect ratio used by 
previous studies is too small for the occurrence of STC. 
We find from our extensive numerical simulation that 
the convective state just above onset is spatiotemporally 
chaotic, which is evident from the snapshot images of 
the vertical velocity field and from the dynamical be- 
havior of three important global quantities: the viscous 
dissipation energy, the thermal dissipation energy and 
the work done by the buoyancy force. This thus pro- 
vides another example of a direct transition to STC, in 
addition to the Kiippers-Lortz transition , ac driven 
electroconvection and a few others . Our method sug- 
gests that by studying the dynamical behavior of some 
global quantities of systems which exhibit STC, one may 
obtain valuable information about the temporal chaos of 
these systems. We also measure the fractal dimensions 
of the global quantities and find a value of about 1.4. 
In addition, we investigate the nature of the conduction 
to STC transition, as well as certain properties of the 
spatiotemporally chaotic state. Our results for the corre- 
lation length suggest that the transition is second order, 
with a mean field power law behavior. In comparison 
recent experimental results for the Kiippers-Lortz tran- 
sition are not consistent with a mean field power law be- 
havior 01 . We present below a simple but rather accurate 
theory for the behavior of the time-averaged convective 
current as a function of the reduced control parameter. 
Finally, we show that the time-averaged structure fac- 
tor (power spectrum) exhibits a scaling behavior with 
respect to the correlation length similar to that found in 
critical phenomena. 

The Boussinesq equations, which describe the evolu- 
tion of the velocity field u(x,y,z,t) — {u,v,w) and the 
deviation of the temperature field 9{x,y,z,t) from the 
conductive profile, can be written in dimensionless form 
as 

V • w = 0, 

du/dt-\-{u-V)u^~Vp + aee^-\-aV^u, (1) 
d9/dt + u-Ve = \7^e + wR, 

where is the unit vector in the vertical z-direction. 
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In the idealized limit of a laterally infinite system, the 
critical Rayleigh number Rc = 277r^/4 and the onset 
wavenumber kc — 7r/\/2- The efficient Marker-and-Cell 
(MAC) 1 15 1^ numerical technique is employed. (In de- 
veloping our algorithm, we have carried out two tests: 
One is against the theoretical results of Schliiter et al. 

the other is against the numerical results of Kir- 
chartz and Oertel In both cases agreements within 
1% have been found.) The boundary conditions for the 
velocities u are free-slip at the upper and lower surfaces, 
and no-slip at the sidewalls. The temperature deviation 9 
is zero on the top and bottom surfaces, while the temper- 
ature gradient W9 normal to the sidewall is set to zero. 
For the initial conditions of u and 9, we have tried both 
random values and values of Gaussian distribution, inside 
a possible range. Since no difference has been found, the 
actual calculation is carried out with initial conditions 
of Gaussian distribution. Our parameters are cr = 0.5 
and 0.03 < e < 0.5, where e — {R — Rc)/Rc is the re- 
duced Rayleigh number. We use mesh points N^jX NyX 
=256 X 256 xl8 and a grid size Ax=A?/=60/256, 
Az=l/18 for an aspect ratio F =60 in the simulation. We 
have run for 360 vertical diffusion times before collecting 
data. Considering that the normal relaxation time to ap- 
proach a steady state is about 10 vertical diffusion times 
[Treiax = 2(1 -|- cr)/37r^cre] , we believe that we are well 
beyond any transient regime. 




FIG. 1. A typical image of the spatially disorganized pat- 
tern in the cell. Dark regions correspond to hot rising fluid 
and white regions correspond to cold descending fluid. The 
vertical velocity fleld u;(x,y,z=l/2) for e = 0.1 is shown here. 

For the low Prandtl number fluid studied here, the con- 
vective pattern near onset has an irregular space-time de- 
pendence. In Figure 1 we show a snap shot image of the 
vertical velocity field w{x, y, z — 1/2) from the numerical 
simulation at e = 0.1. In this image, the apparently dis- 
organized spatial pattern consists of superimposed rolls 



with many different orientations. The time evolution of 
these rolls is through an interface motion, which main- 
tains the type of spatial disorder shown in Figure 1 . Sim- 
ilar images are found for other values of e. It is obvious 
from such images that the convection pattern near onset 
is random in space. 

To illustrate the temporal chaos of the system, we now 
investigate the dynamics of the global quantities which 
characterize the underlying physics of Rayleigh-Benard 
convection. Using (/) to denote the average of / over 
the whole system and taking into account the boundary 
conditions as well as the incompressibility condition, we 
obtain ^d{u ■ u)/dt = F2{t) - Fi{t), and \d{9'^)/dt = 
F^it) - Fait), where (a) Fi = ^cr{{dui/dxj + duj/dxif) 
is the kinetic energy dissipated by the viscosity; (b) 
F2 — (j{w9) is the work done by the buoyancy force; (c) 
F3 = (V6' • V^) is the dissipative thermal energy (gen- 
eration of entropy) owing to temperature fluctuations; 
and (d) F4 = R{w9) = (i^Fj/cr) is the flow of the en- 
tropy fluctuations carried by the vertical velocity. It is 
clear that in the special case of steady state (^ = 0), 
one recovers the condition Fi — F2 and ^4 = ^3. These 
global quantities provide us with a complete description 
of "energy-balance" in the Rayleigh-Benard system. 



0.2 



Fi(t) 



0.19 
0.18 
0.17 
0.16 
0.15 
0.14 
0.13 
0.12 



Fi- 
F2- 
F3- 



Ml 



1 1 



50 



100 150 200 250 
time 



350 400 450 



FIG. 2. A plot of global quantities Fi{t), F2{t) and F3{t) 
as functions of time for e = 0.2. Note that they lie on top 
of each other with differences, though substantial, too small 
to be seen. The time is in units of vertical thermal diffusion 
time ty — d^/K and the origin corresponds to t = 3S0ty in 
real calculation. 

We plot a representative time series of these quanti- 
ties, Fi(t), F2it) and Fait), in Figure 2 for e = 0.2. 
We have rescaled Fi, F2 by aRc, and F3 by RRc so 
that we have Fi = F2 = F3 in a steady state. (Note 
that F4 is simply related to F2 by the factor R/cr.) The 
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most important implication of this figure is the apparent 
chaotic behavior of these quantities over the time interval 
that is accessible to us. To be more concrete, we apply 
the Grassberger and Procaccia method to compute 
the fractal dimensions Df for these quantities and find 
Df ^ 1.42 ± 0.02. This of course is different from the 
fractal dimension that is normally used to characterize 
STC, which diverges with the system size. We believe 
that such global quantities might provide a relatively sim- 
ple way to characterize the temporally chaotic nature of 
spatiotemporally chaotic states such as studied here, al- 
though data over a longer time interval will be necessary 
for such an analysis. It is also interesting to observe that 
the dynamics of these three quantities are almost exactly 
the same, i.e. Fi{t) ~ F2{t) ~ ^3(i), as shown in Figure 
2. Although the differences among them are substantial 
and beyond numerical uncertainties, they are too small 
to be seen under the scale of Figure 2. In fact, this is 
the case for all e studied in the range 0.03 < e < 0.5. 
This is certainly a surprising result considering the ir- 
regular spatiotemporal state we observed. However, a 
theoretical understanding of this has been obtained, as 
outlined later. This result also implies that the quantity 
= cri^a/ (2F2 — Fi), which is often used as a variational 
formulation to determine the critical Rayleigh number, 
behaves as if the system is almost in a steady state. 
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FIG. 3. A plot of kS{k)/(^ vs X — {k ~ kmax)^ (in units of 
fcc), showing scaling and the scaling function F(x) defined in 
the text. Insert: The time-averaged function kS{k) vs. k/kc 
for e = 0.03, 0.05 and 0.1. 

In order to gain more insight into the nature of the 
transition to STC near onset, we have studied the two- 
dimensional structure factor (Fourier power spectrum). 
Since the snapshot images of the patterns appears to 



be isotropic azimuthally, we calculate the azimuthally 
averaged structure factor, and then average the images 
over time, to obtain the time-averaged structure factor 
S{k). The function kS{k) is shown for several different 
values of e in the insert in Figure 3. We also show in 
this figure that kS{k) satisfies a scaling behavior some- 
what similar to that found in critical phenomena, namely, 
kS{k)/£^ = F[{k — kmax)S,], where ^ is the correlation 
length (defined below) and where we have normalized the 
integral of S{k) over fc-space to be unity. Here k^ax is the 
wavenumber which maximizes kS{k), which we choose to 
give the best fit to scaling. Another interesting feature 
of the structure factor is associated with the power-law 
behavior of S(k) ~ k^^ for large wavenumber. This fea- 
ture is observed over the range of e studied here. It is 
interesting to note that this is the same power-law behav- 
ior observed in phase separating systems, in two dimen- 
sions, where it is known as Porod's law. In both cases it 
results from the linear behavior of the real space corre- 
lation function, C(r), for small r, where this correlation 
function is the azimuthal average of the inverse Fourier 
transform of the structure factor. 




FIG. 4. A plot of ^ vs. e. The vertical bars indicate 
the standard deviation, and the solid line corresponds to 
C'^ = ^(7^6 with ^0 = 0.78. These data can also be fitted 
with a nonconventional exponent: see the text for more de- 
tail. 



We have also calculated the correlation length ^ as a 
function of the control parameter e, where we define the 
correlation length through the variance of the wavenum- 
ber, i.e. ^ = (/c^ ~ k )^^^^. The moment /c" is defined 
as F = / \k\'^S{k)d^k/ J S(k)(Pk and S{k) is the time- 
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averaged structure factor. We find that the correlation 
length ^ appears to diverge as e approaches the transi- 
tion point, with a power-law behavior of ^ = Co(e — ^c)^'' 
with z/ = 0.472 ± 0.016, Co = 0.82 ± 0.04 and = 0.005. 
(The fact that Cc is finite instead of zero is due to finite 
size effects.) The behavior of the correlation length is 
also consistent with a mean field power low exponent of 

= 0.5 and = 0.78. This is illustrated in Figure 4. 
The amplitude value is a factor of 3/2 larger than the 
value Co = \/8/37r2 = 0.52 calculated from the curvature 
of the marginal stability curve. 

In order to investigate the heat transport in STC 
near onset, we calculate the Nusselt number N{t) — 
1 + {wff) /i?, which describes the ratio between the heat 
transport with and without convection, as a function of e. 
The quantity A'^ — 1 can be fit with a power law behavior 
of the form iV - 1 = (e - t^Y /g with /i = 1.034 ± 0.025, 
Cc = 0.012 and g = 1.27 ± 0.03. Figure 5 shows that the 
time-averaged — 1 is also consistent with a linear re- 
lation (e — ec)/5, where ec — 0.01 and g — 1.23. (Again, 
owing to finite size effects, the value of ec is nonzero.) We 
have also determined the time-averaged vertical "vortex 
energy" 17 = (w^), where lUz is the vertical vorticity, as a 
function of e from our simulations and observed a power- 
law behavior of 51 ^ with A = 5/2. 
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FIG. 5. Time- averaged TV — 1 vs e near onset. The vertical 
bars indicate the standard deviation, and the solid line is the 
fit of (e — ^c)/'g to the data with ec = 0.01 and g = 1.23. A 
slightly diflferent fitting form is given in the text. 

For theoretical understanding of some of our results, 
we notice that the velocity u and the temperature deriva- 
tion 9 near onset can be approximated by an order pa- 
rameter 'ip{r) in two-dimensional space f multiplied by 
known prefactors with z-dependence pO| . It is then 
straightforward to rewrite the global quantities Fi, F2 



and and the Nusselt number N in terms oi '4){r). We 
further assume that only those modes inside the vicinity 
of kc are excited and equally excited. From these approx- 
imations, we confirm (after rescaling mentioned earlier) 
that Fi{t) ~ F2{t) ~ Fz{t). We also obtain that 

9 = 5(-l) + (2/7^) / 3(cos a)da 
Jo 

= 0.855951 + 0.0458145(7"^ + 0.0709325f7-^ (2) 

where a is the angle between k and some reference di- 
rection, and we have used the explicit form of g(cosa) 
given by Schliiter et al. ||l7|] for free-free boundary con- 
ditions ||2^. For (7=0.5, we find g= 1.2313, which is in 
surprisingly good agreement with the numerical results. 
The theory of the vortex energy ~ e^/^ is more com- 
plicated than the above. All this theoretical analysis will 
be presented elsewhere. 

In summary, we have presented a large scale numeri- 
cal simulation of pattern formation in three dimensional 
Rayleigh-Benard convection. We have calculated the spa- 
tial correlation length and the Nusselt number as a func- 
tion of the reduced control parameter, as well as the dy- 
namics of the viscous energy, the dissipative thermal en- 
ergy and the work done by the buoyancy force. Our 
numerical studies suggest that the transition from the 
conduction state to spatiotemporal chaotic state near on- 
set is a continuous (second order) transition, with mean 
field exponents for the correlation length and the time- 
averaged convective current. We have also demonstrated 
that the time-averaged structure factor satisfies a scaling 
behavior with respect to the correlation length. We be- 
lieve that more studies of scaling in STC by systematic 
experiments for smaller e, larger aspect ratio, and with 
both free- free and rigid-rigid boundary conditions will be 
challenging and valuable. 
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